function lu_pivoting()
% LU_PIVOTING 带选主元的LU分解演示
% 
% 演示部分选主元和完全选主元的LU分解

fprintf('=== 带选主元的LU分解 ===\n\n');

%% 1. 部分选主元的必要性
fprintf('1. 部分选主元的必要性\n');
fprintf('当主元很小时，不选主元的LU分解数值不稳定\n\n');

% 创建小主元矩阵
epsilon = 1e-10;
A = [epsilon, 1, 0; 
     1, 1, 1; 
     0, 1, 1];
b = [1; 3; 2];

fprintf('小主元矩阵 (ε = %e):\n', epsilon);
disp(A);
fprintf('条件数: %.2e\n', cond(A));

% 不选主元的LU分解
try
    [L_no_pivot, U_no_pivot] = lu_decomp_basic(A);
    fprintf('不选主元的LU分解:\n');
    fprintf('L =\n'); disp(L_no_pivot);
    fprintf('U =\n'); disp(U_no_pivot);
    
    % 求解
    y = forward_substitution(L_no_pivot, b);
    x_no_pivot = backward_substitution(U_no_pivot, y);
    fprintf('不选主元的解: x = [%.6e; %.6f; %.6f]\n', x_no_pivot);
    
    residual_no_pivot = A * x_no_pivot - b;
    fprintf('残差: ||A*x - b|| = %e\n', norm(residual_no_pivot));
catch ME
    fprintf('不选主元的LU分解失败: %s\n', ME.message);
end

%% 2. 部分选主元的LU分解
fprintf('\n2. 部分选主元的LU分解\n');
fprintf('PA = LU，其中P是置换矩阵\n\n');

[L_pivot, U_pivot, P] = lu_partial_pivot(A);

fprintf('置换矩阵 P =\n');
disp(P);
fprintf('L矩阵 =\n');
disp(L_pivot);
fprintf('U矩阵 =\n');
disp(U_pivot);

% 验证分解
PA = P * A;
LU_pivot = L_pivot * U_pivot;
fprintf('验证 PA - LU 的误差: %e\n', norm(PA - LU_pivot, 'fro'));

% 求解 PA*x = Pb
Pb = P * b;
y_pivot = forward_substitution(L_pivot, Pb);
x_pivot = backward_substitution(U_pivot, y_pivot);
fprintf('选主元的解: x = [%.6f; %.6f; %.6f]\n', x_pivot);

residual_pivot = A * x_pivot - b;
fprintf('残差: ||A*x - b|| = %e\n', norm(residual_pivot));

% 与MATLAB内置函数比较
[L_matlab, U_matlab, P_matlab] = lu(A);
x_matlab = U_matlab \ (L_matlab \ (P_matlab * b));
fprintf('MATLAB解: x = [%.6f; %.6f; %.6f]\n', x_matlab);
fprintf('差异: ||x_pivot - x_matlab|| = %e\n', norm(x_pivot - x_matlab));

%% 3. 增长因子分析
fprintf('\n3. 增长因子分析\n');
fprintf('增长因子衡量LU分解过程中元素的增长\n\n');

% 计算增长因子
max_A = max(abs(A(:)));
max_U = max(abs(U_pivot(:)));
growth_factor = max_U / max_A;

fprintf('原矩阵最大元素: %.6f\n', max_A);
fprintf('U矩阵最大元素: %.6f\n', max_U);
fprintf('增长因子: %.6f\n', growth_factor);

% 最坏情况的增长因子示例
fprintf('\n最坏情况示例 (增长因子 = 2^(n-1)):\n');
n = 5;
A_worst = eye(n) + tril(ones(n, n), -1);
A_worst(:, end) = 1;
A_worst(end, :) = [-1, -1, -1, -1, 1];

fprintf('最坏情况矩阵 (%d×%d):\n', n, n);
disp(A_worst);

[L_worst, U_worst, P_worst] = lu_partial_pivot(A_worst);
max_A_worst = max(abs(A_worst(:)));
max_U_worst = max(abs(U_worst(:)));
growth_worst = max_U_worst / max_A_worst;

fprintf('增长因子: %.1f (理论最大值: %.1f)\n', growth_worst, 2^(n-1));

%% 4. 完全选主元
fprintf('\n4. 完全选主元\n');
fprintf('PAQ = LU，其中P和Q都是置换矩阵\n\n');

% 创建测试矩阵
A_complete = [1e-10, 1, 2; 
              2, 1e-10, 1; 
              1, 2, 1e-10];

fprintf('测试矩阵:\n');
disp(A_complete);

[L_complete, U_complete, P_complete, Q_complete] = lu_complete_pivot(A_complete);

fprintf('行置换矩阵 P =\n');
disp(P_complete);
fprintf('列置换矩阵 Q =\n');
disp(Q_complete);
fprintf('L矩阵 =\n');
disp(L_complete);
fprintf('U矩阵 =\n');
disp(U_complete);

% 验证分解
PAQ = P_complete * A_complete * Q_complete;
LU_complete = L_complete * U_complete;
fprintf('验证 PAQ - LU 的误差: %e\n', norm(PAQ - LU_complete, 'fro'));

%% 5. 性能比较
fprintf('\n5. 性能比较\n');
fprintf('比较不同选主元策略的性能\n\n');

sizes = [100, 200, 400];
fprintf('矩阵规模    无选主元    部分选主元    完全选主元\n');
fprintf('--------    --------    ----------    ----------\n');

for i = 1:length(sizes)
    n = sizes(i);
    A_test = randn(n, n) + eye(n);
    
    % 无选主元
    tic;
    try
        [L1, U1] = lu_decomp_basic(A_test);
        time_no_pivot = toc;
    catch
        time_no_pivot = inf;
    end
    
    % 部分选主元
    tic;
    [L2, U2, P2] = lu_partial_pivot(A_test);
    time_partial = toc;
    
    % 完全选主元
    tic;
    [L3, U3, P3, Q3] = lu_complete_pivot(A_test);
    time_complete = toc;
    
    fprintf('%6d      %6.3f s      %6.3f s      %6.3f s\n', ...
            n, time_no_pivot, time_partial, time_complete);
end

%% 6. 稳定性测试
fprintf('\n6. 稳定性测试\n');
fprintf('使用随机病态矩阵测试数值稳定性\n\n');

% 生成病态矩阵
n = 10;
U_test = triu(randn(n, n));
L_test = tril(randn(n, n));
L_test = L_test - diag(diag(L_test)) + eye(n);  % 单位下三角

% 构造病态矩阵
D = diag(logspace(-10, 0, n));  % 条件数很大
A_ill = L_test * D * U_test;
b_ill = randn(n, 1);

fprintf('病态矩阵条件数: %.2e\n', cond(A_ill));

% 部分选主元求解
[L_ill, U_ill, P_ill] = lu_partial_pivot(A_ill);
Pb_ill = P_ill * b_ill;
y_ill = forward_substitution(L_ill, Pb_ill);
x_ill = backward_substitution(U_ill, y_ill);

% 计算残差和误差
residual_ill = A_ill * x_ill - b_ill;
x_exact = A_ill \ b_ill;
error_ill = norm(x_ill - x_exact) / norm(x_exact);

fprintf('相对残差: %e\n', norm(residual_ill) / (norm(A_ill) * norm(x_ill)));
fprintf('相对误差: %e\n', error_ill);

end

function [L, U] = lu_decomp_basic(A)
% 基本LU分解（无选主元）
n = size(A, 1);
L = eye(n);
U = A;

for k = 1:n-1
    if abs(U(k, k)) < eps
        error('主元为零');
    end
    
    for i = k+1:n
        L(i, k) = U(i, k) / U(k, k);
        U(i, k:n) = U(i, k:n) - L(i, k) * U(k, k:n);
    end
end
end

function [L, U, P] = lu_partial_pivot(A)
% 部分选主元的LU分解
n = size(A, 1);
L = eye(n);
U = A;
P = eye(n);

for k = 1:n-1
    % 找到第k列中绝对值最大的元素
    [~, max_idx] = max(abs(U(k:n, k)));
    max_idx = max_idx + k - 1;
    
    % 交换行
    if max_idx ~= k
        U([k, max_idx], :) = U([max_idx, k], :);
        P([k, max_idx], :) = P([max_idx, k], :);
        if k > 1
            L([k, max_idx], 1:k-1) = L([max_idx, k], 1:k-1);
        end
    end
    
    % 消元
    for i = k+1:n
        L(i, k) = U(i, k) / U(k, k);
        U(i, k:n) = U(i, k:n) - L(i, k) * U(k, k:n);
    end
end
end

function [L, U, P, Q] = lu_complete_pivot(A)
% 完全选主元的LU分解
n = size(A, 1);
L = eye(n);
U = A;
P = eye(n);
Q = eye(n);

for k = 1:n-1
    % 找到剩余子矩阵中绝对值最大的元素
    [max_val, max_idx] = max(abs(U(k:n, k:n)), [], 'all', 'linear');
    [row_idx, col_idx] = ind2sub([n-k+1, n-k+1], max_idx);
    row_idx = row_idx + k - 1;
    col_idx = col_idx + k - 1;
    
    % 交换行
    if row_idx ~= k
        U([k, row_idx], :) = U([row_idx, k], :);
        P([k, row_idx], :) = P([row_idx, k], :);
        if k > 1
            L([k, row_idx], 1:k-1) = L([row_idx, k], 1:k-1);
        end
    end
    
    % 交换列
    if col_idx ~= k
        U(:, [k, col_idx]) = U(:, [col_idx, k]);
        Q(:, [k, col_idx]) = Q(:, [col_idx, k]);
    end
    
    % 消元
    for i = k+1:n
        L(i, k) = U(i, k) / U(k, k);
        U(i, k:n) = U(i, k:n) - L(i, k) * U(k, k:n);
    end
end
end

function x = forward_substitution(L, b)
% 前向替换
n = size(L, 1);
x = zeros(n, 1);
for i = 1:n
    x(i) = (b(i) - L(i, 1:i-1) * x(1:i-1)) / L(i, i);
end
end

function x = backward_substitution(U, b)
% 后向替换
n = size(U, 1);
x = zeros(n, 1);
for i = n:-1:1
    x(i) = (b(i) - U(i, i+1:n) * x(i+1:n)) / U(i, i);
end
end